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ABSTRACT 

We consider weak solutions of hyperbolic systems in primitive (non-conservation) form for 
which a consistent conservation form exists. We show that for primitive formulations, shock 
relations are not uniquely defined by the states to either side of the shock but also depend 
on the viscous path connecting the two. Scheme- dependent high order correction terms are 
derived that enforce consistent viscous shock profiles. The resulting primitive algorithm is 
conservative to the order of the approximation. One dimensional Euler calculations of flows 
containing strong shocks clearly show that conservation errors in primitive calculations are 
reduced to truncation levels and that both conservative and primitive flow calculations are 
of comparable quality. 


1 Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1. INTRODUCTION 


It is common wisdom that consistent shocked solutions can be numerically 
captured only if the numerical algorithm meets discrete conservation 
requirements. Indeed, conservative numerical calculations satisfyling the 
entropy condition converge to the correct physical solutions as the mesh 
size tends to zero [4]. Conversely, examples easily show that calculations 
which are not conservative converge to completely non-physical solutions- 
(eg [5]). Non-conservative or primitive formulations, however, are 
strikingly simpler and are less coupled than their conservative 
counterparts. As such they may offer advantages, either in computational 
efficiency or in accuracy gains. In Fluid Dynamics, the diagonal 
characteristic formulation is probably the most prominent example [6] . 
Other formulations using the entropy function are also favoured by many. 
The less common choice of velocity components and pressure, often leads to 
accuracy gains near contact surfaces separating materials of different 
types [3]. In multi-dimensional setups, this choice also enables the 
advection of a uniform passive velocity, completely decoupled from the ID 
Riemann solver in the cross direction. In high speed near vacuum 
conditions, the internal energy is an important but usually very small 
quantity, which due to numerical truncation errors may become negative. 
This problem may be resolved by using internal energy as a dependent 
variable at the cost of sacrificing conservation. Exaples of primitive 
forms arising in Elasticity are discussed by Colombeaux and Le Roux in 
[ 1 , 2 ]. 

Straight forward discretizations of primitive formulations result in both 
incorrect shock speed/location and wrong jump across shock transition. We 
show, that primitive formulations do not possess unique jump conditions for 
steady viscous shock profiles. Jump conditions depend not only on the 
limiting left and right states but also on the viscous path that connects 
them. This has also been shown in [1,2] using arguments from generalised 
functions theory. The secret of correct shock capturing thus lies in 
getting the path right. Although physically, there is only one correct 
such path, numerically there are many. In fact, as many as there are 
conservative numerical schemes. Indeed, while physical shock transition is 
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governed by physical viscosity mechanisms, numerical shock transition is 
governed by numerical viscosity mechanisms, whose precise form depend on 
numerical truncation errors. The analysis of Le Roux and Colombeaux in 
[1,2] tries to enforce physical microscopic behaviour on the numerical 
algorithm. The physical microscopic behaviour is deduced either from a 
consistent conservation form or from empirical data. Ignoring numerical 
viscosity mechanisms, this analysis is not fully justifiable on the 
discrete level. In contrast, the analysis presented in this work is 
performed directly on the dicrete level and enforces correct numerical 
microscopic behaviour. We follow an idea introduced by Zwas and Roseman 
[10] , who looked into the effect of non-linear transformations on weak 
solutions of conservation laws. They have considered the particular case 
of an original set of conservation laws which transforms into another set 
of conservation laws. They have looked at the viscous fora of the 
equations and showed that unless the viscosity terms are included in the 
transformation, the latter set of conservation laws will produce 
inconsistent weak solutions. For a conservative system to transform into 
another pseudo-conservative system is, however, a very special case. More 
commonly, it transforms into a primitive form. . Shocks obtained by 
primitive calculations depend inherently on getting the underlying viscous 
path right. We follow [10] in the case of general formulations written in 
primitive form, for which a consistent underlying conservation form exists. 
This is where we believe its great promise rests. We consider hyperbolic 
primitive formulations 


v + A (w)w = 0 

-t - -X 

and derive general, scheme dependent, high order correction terms 

w + A (w)v = Lt v f(w,w ,w ,\) 

— t - “X — “X — t 

where p is the order of the scheme and X=At/Ax the mesh ratio. Their 
inclusion on the RHS of the primitive formulation renders the viscous forms 
of the conservative and primitive algorithms equivalent. Though not 
strictly conservative, the resulting primitive algorithm is conservative to 
the order of the approximation. Correction terms are obtained for the 
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first order Lax-Friedrichs and upwind schemes without reference to a 
particular system. Their specific form is given for the ID Isothermal Euler 
equations and the complete ID Euler equations. The effect of the 
correction terms is demonstrated on one dimensional Euler calculations of 
flows containing strong shocks. It is clearly seen that errors in weak 
solutions are reduced to truncation levels, and that both conservative and 
primitive flew calculations are of comparable quality. 


2. WEAK SOLUTIONS AND VISCOSITY 

Consider scalar conservation laws described by the Initial Value Problem 
(IVP) , 


u + f(u) = 0 (1) 

l X 

u (x, 0) = (x) 

Denote by a(u) = df/du the characteristic speed of the equation, then 
solutions to (1) can be written implicitly as 

u(x,t) = <p(x-a(u)t) 

Depending on whether a(<p(x)) is an increasing or decreasing function of x, 
an initially smooth solution u(x,t) will either remain smooth or develope 
discontinuities or shocks. Integral conservation considerations allow the 
solution to be extended beyond the time of shock formation. The broader 
concept of Veak Solutions is introduced, describing piecewise smooth 
solutions separated by curves of discontinuity, across which the solution 
satisfies the Rankine-Huganiot jump conditions 


s 



(2) 


Here s is the shock speed and ( ) R and ( ) L denote the states to its 
immediate right and left. Weak solutions, however, are not unique. The 
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criterion that rules out all but one of the solutions is known as the 
Entropy Condition. This condition can be shown [5] to select a unique 
solution which is the limit of solutions u^(x,t) of the Viscous problem 

u + f (u) = eu 

t X XX 

as viscosity vanishes c -* 0. The concept of viscosity, thus lies in the 
heart of correct, entropy satisfying, shock representation. While in 
smooth parts of the flow, the viscosity term can be neglected on grounds of 
order of magnitude, in regions of rapidly varying solutions its neglect 
leads to ambiguities. These can only be resolved upon conceptual 

re-introduction of the neglected terms. 

CONSERVATION FORMS. PRIMITIVE FORMS AND VISCOUS SHOCK PROFILES 


The more general viscous form of the equation reads 


u + f(u) = e(F(u)u ) (3) 

1 X XX 

for some function F(u). Assume that F(u) is such that stable shock 
profiles exists and consider a steady viscous shock profile moving at a 
constant speed s, u = u (x-st), satisfying 


u 


» u 

X-* -<» L 


U 


U — H > 

x X “ » ±00 


x — ► +00 
0 


u 


R 


Substituting (4) into (3) and integrating over xe(-<x>, ou) gives 


(4) 


-s(u - u ) + (f - f ) = l [(F(u)u ') - (F(u)u') 1 (5) 

K L h. L K L 

By (4c), the RHS of (5) vanishes, yielding the jump conditions (2). 
Provided F(u) is admissible, this result does not depend on its precise 
definition although the viscous path connecting u and u obviously does. 

Let a transformation T be defined by Tdu=dw and assume there exists a g (v) 
such that Tdf (u) =dg (w) . Then (1) transforms into the conservation law. 
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( 6 ) 


w + g(w) ~ 0 

1 X 


with shock solutions satisfying 


s 


* g 


L 


- W 

I. 


which is inconsistent with (2). Any other admissible non-zero RHS in (6) 
of the form t(G(w)w ) yields the same jump relations. If, however, 

X X 

equation (3) is transformed together with the viscosity term in (3), then 
it reads 


v + g(") = eT(F(u)u ) (7) 

t X XX 

= e (TF(u)u ) - CT F(u)u 

XX X X 

Substituting a viscous shocK profile (4) into (7) now yields 

OD 

s(ff„ - v .) = (g - a ) + f T F(u)u dx (8) 

K L K L X X 

-CD 

and should give correct shock speed provided the viscous profile used in 
the integration is consistent with (3). 

Note that the transformed equation (6) may not always be written in 
conservation form. More generally it reads 

w + b(w)« = 0 (9) 

t X 

for some b(w ). Note that b(w) always satisfies the conservation law 

i>(w) t + [\ (b(w)) 2 ] = 0 (10) 

1 Z x 

but this in itself does not make (10) more correct than (6) . That will 
depend on whether the assumed underlying viscous form is the correct one. 
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Viscous conservative systems read 


u + f(u) = t{F(u)u ) 

~t - - X ““XX 

where F(u) is now a Viscosity Matrix. Again, not every matrix F(u) yields 
stable shock profiles. Under certain assumptions, the identity matrix I is 
admissible (see for example [7] and references cited therein) . 

Primitive formulations of hyperbolic systems depend crucially on the 
correct choice of viscous paths. Primitive systems have the form 

w + A (w)w - 0 

- t - — X 

where A(w) is not a Jacobian matrix with respect to w. If the viscous 
primitive form is assumed to be 

w + A(w)w = cw ( 11 ) 

— t — “X “XX 

steady viscous shock profiles satisfy 

w 

— R 

/ A(ff)dw = s(w^-w^) (12) 

w 

-L 

Since A (w) is not a Jacobian with respect to w, the integration is path 
dependent and so are both the shock speed s and the jump (w - w ) (see 
Figure (1)). They will only be correct if the integration is along a 
consistent path, ie if a consistent RHS is taken in (11). If w varies 
linearly across shock transition, then by (12) 

x 

R 

( X 1 — / = S(W -*~ W -J 

R L X L 

implying that the jump (v -w ) is an eigenvector of a path-dependent 
average of A(v) , and that the shock speed s is the associated eigenvalue . 
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4. NUMERICAL SOLUTIONS, NON-LINEAR TRANSFORMATIONS AND CORRECT SHOCK 


REPRESENTATION 


Attempting to solve either (1) or (9) numerically immediately raises the 
question of consistent viscous integration paths since due to numerical 
viscosity, captured shocks always get smeared over a number of grid points.- 
The precise form of shock transition depends on numerical truncation 
errors. Consequently, their relevance to physics is not in their precise 
details but in some average interpretation of shock location and in an 
asymptotic interpretation of the limiting states to either of its sides. 
While there is only one correct physical shock transition, there are many 
correct numerical shock transitions. Indeed, let a numerical grid be 
defined by the partition parameters (Ax, At) and let u" as u(jAx,nAt). 
Consider the system of conservation laws. 


u + f(u ) = 0 (13) 

-t - - X 


Then any conservative numerical scheme 


n + 1 n 

a ~ U 
~j “j 




that consistently approximates (13) produces shock transitions which are 
correct in that average sense. Here h" = h(u p , — ,u . ) is the 

-j + l/2 - -j --C+ t “ j ♦ / l 

numerical flux function at the j + i/2 cell interface, with £ and a denoting 
the numerical stencil width, and consistency implies 


h n (u,u , . . • ,u) = f(u) 

-j + 1 / 2 - - ” -- 

Other considerations dictate which shock representation is more acceptable. 
The role of the viscous path is revealed in a more concrete way by writing 
the viscous form associated with a given numerical scheme [10] . Keeping the 
leading order terms in the numerical truncation error this reads 


(14) 


n+ 1 
U 

~j 



+ U J 

-J + l 


5 «” , 

2 -J + l 


C.' 


with the numerical flux function 


h = 3 ff + O 

-j + 1 / 2 Z “j + 1 -j 


1 

21 


/ n 
(U . 

“J + 


1 



the viscous form reads 

u + f(u) = (u / \ Z - u ) (15) 

and since (14) is conservative, the resulting viscous path is consistent, 
though not unique. Let w be a different set of dependent variables, let T 
the Jacobian of the transformation Tu - dw. Premultiplying (13) by T 
yields the primitive system. 


w + A (w)w = 0 

“t — “ X 


(16) 


Let (16) be approximated by a 'LxF- type’ approximation 


n ^ 1 1 , n n , 1 , , . . » , n ^ i 

v. = - (w. + w. ) - (A. + A. )(«. - v. ) 

-j 2 -j-1 -j+l 4 j + l j-i -j + l -j-t 


(17) 


Then its viscous form reads 


w + A ( v) w 

_t - _ x 


At 

2 


(w / X 2 - w ) 

—XX -It 


(18) 


Unless this viscous form is equivalent to (15) , or for that matter to the 
viscous form of any other first order conservative scheme, it will yield 
inconsistent weak solutions. In other words, let 


D = T(u A 2 - u ) 

-XX “ t t 


(w A Z 


~ E ) 


(19) 


then (13) and (16) will converge to the same weak solutions (to order At) 
only if D = 0. This requirement will not in general be met. To enforce 
correct weak solutions, the correction terms D must be added to the RHS of 
the primitive formulation (14), which should read instead. 
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(19) 


w + A ( «) w 

-t — — x 


At 

2 


D 


5. EXAMPLES 


The ID isothermal Euler equations are given by 


/ > 
l ) 


f pu 


+ 

z z 


1 

L pu +pc J 


= 0 


( 20 ) 


J 

Here p and u are density and velocity and c = 1 is the constant sound 
speed. A right moving shoch with u^ = 0 satisfies the jump relations. 


/ V. 


s * 


( 21 ) 


Multiplying (20) by the transformation matrix 

T = (-U/P 1/p ] 

gives a primitive formulation in terms of p and u 


f > 

p 


f 

u 

p 

P 


+ 

1 



u U ; 


( 22 ) 


Although equation (23) may be rewritten in pseudo-conservation form 


/■ > 
P 


pu 


+ 

1 

u Z /2 +ln(p) J 


(23) 


it does not represent any genuine physical conservation and will give 
non-physical weak solutions. Indeed, the jump relations for system (23) 
read (compare (21)), 
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2 

U 

L 


( 24 ) 


2 ( P-PJ 

= — ^-r/ln (p /p ) 

P L P K } L R 

5 - 


A third formulation, in terms of ln(p) and u takes a symmetric form. 


'in (p)' 

4 

f 

u 

1' 

'In (p/| 

> u J 

t 

1 

u 

> u J 


for which the transformation matrix is 


T = 


[ 1/P 
(“U/p 


° 1 

1/P ) 


(25) 


Denote the above systems by system I, II and III. System I is, in this 
case, the consistent system. The correction terms for systems II and III 
are respectively 


II P 


[p u /X - pul 

V X X 1 


III 


'- i/p rp x P x A - p t p t J] 

J(P«A Z -P«) 

V P XX It 


(26) 


That the first component in Pjj is zero should come as no surprise, since 
this is an equation for the conserved quantity p and requires no correction 
of order At. For computational convenience, the time derivatives in (26) 
may be replaced by spatial derivatives using (22) or (25) . 


The complete ID Euler equations in conservation form read 


f > 

P 


r > 

p v 

pu 

+ 

pu + p 

l E J 

t 

^ uE + up , 


( 27 ) 


Here E is the specific total energy and p the pressure, obtainable from 
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p = (y-D[E - | pu 2 ] 


using the ideal gas assumption. The primitive form using p,u and p reads 


p ] 

u 

+ 

u 

0 

p 

u 

o ' 

i/p 

' p ' 

u 

. p J 


[ 0 

yp 

u j 

. p . 


is obtained by the transformation 


T = 


1 

-u/p 

(y-iju 2 

2 


0 

1/p 

-d-Du 


o 1 

0 

( 1 - 1 ) 


1 28 ) 


and the correction terms are 


D = 


2 . _ 2 , 

-( P u A - P u ) 
p xx it 

[ (y-l)p [u u A z - u u ] J 

V XX t t y 


Consider the first order upwind approximation to (13) 


n ♦ 1 

U 
~ J 


U n - X ( (A c )* ,, 9 (u.- U. J + (A c )~ 

-j j-1/2 -j -J-l j+t/2 -J+1 “J 


(U 


u )) 


(29) 


Here, A° = df/du is the Jacobian matrix, (A ) denote its positive and 
negative parts and (~) indicates locally averaged values. The superscript 
c denotes to conservative formulation. The viscous form of the first order 
upwind is 


a + f(u) = ^ ( (\A° \ u ) / X 

-1 - - X 2 -XX 


V" 


(30) 


Let w be a set of primitive variables and let T - dw/du be the Jacobian of 

~~ * 
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the transformation. The viscous form of the primitive formulation reads 


v + A p (v)w = =£ ((\A p \w ) / \ - w )) 

-t - - -x 2 -XX -tt 

The superscript p denotes primitive formulation. Then 

A C = T~ l A P T 
I A c ! = r _1 U p l T 

The correction terms for the first order upwind are 

D= [ T(\A c \u } - (\A p \w ) 1/X - (Tu -w ) 

V -xx -x x )' -t t -t t 

or after rearrangement, 

D = ( T ( T ') ld p lw ]/x - T ( T X ) w (33) 

V XX/ i t 

For the ID Euler equations in the particular set w = (p,u,p), given in (28) 
the correction terms are 


(31) 


(32) 


D = 


1 ( P x U xV ri/C > » X P X C Z + (<P/C)U X + (1/PC)P X P X )C 4 , 

2p{ X — - 4 P t » J 


(34) 


(1-1) ( p W C 3 + (1/c) W< „ , 

~ l x ' 2 P Vt J 


where 


c = I u-cl + 21 ul +1 u^-cl 
C 2 = lu-c! - 21 ul +lu*cl 
c 3 ~ I u-cl + I u+cl 
c^ = I u-cl - I u+cl 
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6. NUMERICAL EXPERIMENTS 


In all the following Figures, the dashed profiles were obtained by a 
conservative calculation, hence consistent solutions. The solid profiles 
were obtained by a primitive calculation. Figures (2) and (3) describe 
experiments with the ID isothermal Euler equations, given in genuine 
conservation form in equation (20) and in two alternative primitive forms 
in equations (22) and (25) . The conservative form was approximated by the 
LxF scheme (14) and the primitive forms by the 'LxF-Type' scheme (17). The 
correction terms are given by equation (26), where time derivatives are 
replaced by spatial derivatives, nodal values are replaced by central 
averages and x derivatives by centered differences. Initial data for this 
test were (P L < U L ) = (0.4, 1.0) and (p r ,u r 7 = (0.1, 0.0) . The data were 

chosen to yield distinctively different jump conditions for the first two 
systems, given by equations (21) and (24) respectively. As is clear from 
the Figures, adding the correction terms to the primitive formulations 
reduces the errors to truncation level. It may also be noticed that the 
error in Figure (3) is slightly larger than that in Figure (2). This may 
be attributed to the fact that the variable ln(p) is very ^sensitive to 
small changes in p in the density range over which the test was conducted 
and that consequently the respective formulation suffers larger truncation 
errors. In Figure (4), the ID Euler, equations, given in conservation form 
by (27) and in primitive form by (28) are approximated by the first order 
upwind scheme (29), using Roe's averages [8] for the conservation form and 
simple arithmetic averages for the primitive form. The correction terms 
are given by equation (34), where again, time derivatives are replaced by 
spatial derivatives, local values are centrally averaged and x derivatives 
are replaced by centered differences. The solution was found not to be 
sensitive to the manner in which the correction terms were approximated. 
Figure (4) depicts Sod's shock-tube problem, with initial data (p , u ,p ) = 
(1.0, 0.0, 1.0) and = (0.125,0.0,0.1) . Again, the correction 

terms reduce the errors to truncation level. Inspection of the correction 
terms in (34) reveals that all the products that appear in them contain 
either u or p or both. Both these derivatives vanish near contact 

X X 

surfaces, indicating that the correction terms only act away from these 
regions. Applying the correction terms cannot thus affect the resolution 
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of contact surfaces. This was exploited in the tests shown in Figures (5) 
and (6) , where the correction terms (34) are used in conjunction with the 
second order upwind scheme and superbee flux limiter [9] . Figure (5) 
depicts Sod's shock-tube problem. Figure (6) depicts a more severe 
shock-tube test, with initial data (*\' u l'^l> = (1.0, 0.0, 1.0) and 

(P„rU a ,pJ = (0.125,0.0,0.1) , leading to a shock wave of pressure ratio 
4:1. Indeed, the crisp representation of the contact surface is not 
damaged in any way by the correction terms while the errors due to shock 
formation are again removed. This novel feature is peculiar to choices of 
primitive formulations that include u and p, both of which are constant 
across contacts. It cannot, in general, be expected of other primitive 
forms . 


7. CONCLUSIONS 


It has been shown that primitive formulations of conservation laws do not 
possess uniquely defined weak solutions. Jump relations across shocks were 
shown to depend not only on the limiting left and right states^ but also on 
the viscous path connecting the two. A technique has been described to 
enforce consistent weak solutions on primitive formulations. The method is 
based on deriving high order correction terms, that render the viscous form 
of the conservative and primitive formulations equivalent. The resulting 
primitive algorithm is conservative to the order of the approximation. The 
explicit form of the correction terms is scheme-dependent. Expressions 
were obtained for the first order LxF and upwind schemes. This technique 
was implemented to the ID Euler equations in problems containing fairly 
strong shocks. It has been demonstrated that applying the correction terms 
reduced conservation errors to truncation levels and that conservative and 
primitive flow calculations were of comparable quality. This method shows 
great promise with other primitive formulations. 


fr 
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Figure (2) - ID Isothermal Euler Equations: 

Dashed line by conservation form (20) 
Solid Line by primitive form (22) 

(A) Without and (B) with correction terms 







Figure (3) - ID Isothermal Euler Equations: 

Dashed line by conservation form (20) 
Solid Line by primitive form (25) 

(A) Without and <B) with correction terms 
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Figure (4) - ID Euler Equations - Sod’s shock tube problem by first order 
upwind scheme: 

Dashed line by conservation form (27) 

Solid Line by primitive form (28) 

(A) Without and (B) with correction terms. 












Figure (6) - ID Euler Equations - Strong shock tube problem by second order 
upwind scheme: 

Dashed line by conservation form (27) 

Solid Line by primitive form (28) 

(A) Without and <B ) with correction terms. 
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